Structural defects in ion crystals by quenching the external potential: the 
inhomogeneous Kibble-Zurek mechanism 
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The non-equilibrium dynamics of an ion chain in a highly anisotropic trap is studied when the 
transverse trap frequency is quenched across the value at which the chain undergoes a continuous 
phase transition from a linear to a zigzag structure. Within Landau theory, an equation for the 
order parameter, corresponding to the transverse size of the zigzag structure, is determined when 
the vibrational motion is damped via laser cooling. The number of structural defects produced 
during a linear quench of the transverse trapping frequency is predicted and verified numerically. It 
is shown to obey the scaling predicted by the Kibble-Zurek mechanism, when extended to take into 
account the spatial inhomogeneities of the ion chain in a linear Paul trap. 
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The non-equilibrium statistical mechanics of long- 
range interacting systems is one of the challenging prob- 
lem in statistical physics Nevertheless, close to a 
continuous phase transition it is sometimes possible to 
use concepts of equilibrium statistical mechanics in or- 
der to make some predictions for the system when the 
value of a control parameter is quenched through the 
critical value. The Kibble-Zurek mechanism (KZM) has 
become a useful paradigm in this arena, accounting for 
a variety of phenomena ranging from the formation of 
massive particles in the early universe Q to the vortex 
formation in superfluid Helium The model applies 
to systems with a continuous phase transitions which is 
well described within Landau theory, and allows one to 
estimate the density of defects which are formed when 
quenching the control field v t across the critical value. 

In a nutshell, by comparing the characteristic time tq 
of change of the control field v t with the relaxation time 
r(V t ) of the system at equilibrium one identifies the 
corresponding freeze-out time scale, t, which separates 
the regime in which the system follows adiabatically the 
quench, from the regime in which the system behaves as 
if the dynamics was frozen out. The correlation length 
£ at Vt(i) then gives the characteristic length over which 
the system remains correlated, and hence the density of 
defects. The KZM prediction of the density of defects 
has been verified numerically and experimentally 

in a variety of systems [7[. Recently the model has been 
extended to describe the quench dynamics in a quantum 
phase transition The standard model of topological 
defects for homogeneous phase transitions breaks down 
whenever the quench is local or the critical control pa- 
rameter and the resulting transition become spatially de- 
pendent 0, flo| . This situation is arguable ubiquitous in 
nature and the main focus of this work. In this Letter 
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FIG. 1: (a) Charge distribution at equilibrium of an ion chain 
in a linear Paul trap for decreasing transverse trapping fre- 
quency (from top to bottom). Due to harmonic trapping, 
the density of ions is larger in the centre, where the repul- 
sion is larger and the zigzag instability is first evident, (b) 
Charge distribution after a quench through the critical trans- 
verse frequency exhibiting both types of structural defects 
(solid lines): kink (light) and antikink (dark). 



we study the out-of-cquilibrium dynamics in a inhomoge- 
neous laser-cooled Wigner crystal, quenched through its 
critical point. Here, we propose a novel system to test the 
KZM in a inhomogeneous phase transition which enjoys 
unprecedent level of experimental control and amenabil- 
ity for defect detection. The Wigner crystal is composed 
of single-charged ions, which are confined in Paul or Pen- 
ning traps and mutually repel via unscreened Coulomb 
In highly anisotropic traps the ions can 
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form a linear chain, which has a mechanical instability to 
a degenerate chain, with a zigzag structure, controlle d by 
the density or by the transverse trap frequency [l2l Il3| . 
The form of the corresponding distribution of charges at 
equilibrium is shown in Fig. HJa), for the case in which 
the ions are confined by a linear Paul trap. At equi- 
librium such instability is a second order phase transi- 
tion (Hill EH. We determine the scaling of the density 
of defects after crossing such transition as a function of 
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the cooling rate and of the quenching rate of the trans- 
verse trap frequency. In particular, space like separated 
regions may develop zigzag structures with different ori- 
entations [16| , some of which are displayed in Fig. QJb) . 
These regions are the analogs of magnetic domains in 
a ferromagnetic material and the interface between do- 
mains is a structural defect. 

The system we consider is composed by N ions of mass 
to, charge Q and coordinates r n — {x n ,y n , z n ) that are 
confined along the cc-axis by a strongly anisotropic, ra- 
dial trap. The Lagrangian describing the dynamics of 
the ions is L = T — V where the kinetic and potential 
energies take the form, respectively, T — ^rn^n^n an d 

v r h m Enlv 2 xl+v?(yl+4)]+E n ^>Q 2 /mr n -r' n \), 

with v and vt the frequency of the axial and the trans- 
verse confinement. At sufficiently low temperature and 
sufficiently large values of vt the ions form a chain along 
the x-axis. In the local density approximation the lin- 
ear density n(x) is approximated by the function n(x) = 

(l — j^j with L the half-length of the chain and x 

the distance from the center [13] • The interparticle spac- 
ing given by a{x) — l/n(x) is a slowly- varying function 
of the position. In the thermodynamic limit, in which 
a(0) — > a as the number of particles N — > oo, one recov- 
ers in the center of the trap the statistical mechanics and 
dynamical properties of an infinite chain with uniform 
interparticle distance a 12, ljj . The mechanical stabil- 
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ity of the ion chain is warranted provided that v t > u \ , 



such that the critical value 



is a function of the ax- 
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ial trap frequency and the ion density. At v\ the chain 
undergoes a transition to a zigzag configuration, with 
transverse size b [l!| . 

In the thermodynamic limit, the structural change is 
a second order phase transition, with v t a control field 
and b the order parameter 12J. In particular, the critical 
value of the transverse frequency in the thermodynamic 
limit is given by v[ c) = w \AC(3)/2 = (2.051 . . . )uj q , C 
being the Riemann-zeta function, ujq — \J Q 2 /ma 3 and 
a the uniform interparticle spacing. For N finite, the 
transition from a linear to a zigzag chain occurs at the 
value ps 3Nv/ (4-y/log N), where the corrections scale 
with powers of 1/ log N [18j , and corresponds to the rela- 
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s t — i ma ^ 3 with a(0) = l/n(0). For the system 
under harmonic axial confinement, the zigzag instability 
occurs first at the center of the trap, where the density 
is larger, and at lower values of v t extends towards the 
edge of the chain, as sketched in Fig. [Ha). 

We now determine a Ginzburg-Landau (GL) equa- 
tion for the position-dependent transverse size b(x) of 
the zigzag chain, which is here a continuum field mode 
^(x). The derivation is obtained assuming that a coarse- 
grained length scale Sx can be defined, with Sx 3> a(x) 
and a(x) S> \Sx(da(x)/dx)\. In this limit, we make 
a slowly-varying ansatz for the short-wavelength eigen- 



modes of the linear chain, such that for a given eigen- 
mode we can write a n = a n e lkna , with a n slowly- varying 
amplitude 18]. Within the local-density approximation, 
we identify a local value of the critical transverse fre- 
quency v[ c \x) 2 — ^ m a( x yi an d write the Lagrangian 
L = j dxC(x) with the Lagrangian density 
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where ^p a (x) is the field, giving the zigzag size as a 
function of the position, p(x) = mn(x) is the linear 
mass density, and S(x) — vf — v\ (x). The parame- 
ter h(x) = ojQa(x) Vl°g 2 is a velocity, and determines 
the speed with which a transverse perturbation propa- 
gates through the chain. Finally, the parameter A(x) = 
(93C(5)/32)wQ/a(a;) 2 is positive and determines the value 
of the order parameter where 5{x) < 0. The Lagrangian 
density we derived has the form of a GL equation. It is 
valid for the modes of the linear chain close to instability 
and extends the theory presented in [l2T ]. 

The minimal energy solution of Eq. (p} fulfills the re- 
lation: ^(Six) + 2A(x)[(ip y ) 2 + {V) 2 ]) = 0, which al- 
ways admits the solution tfj a — corresponding to all 
the ions on the x axis, and it is stable only for 5(x) > 
for all x, i.e. in the linear chain phase. For 5(x) < 
there is a continuous manifold of solutions of the form, 
g(x) = y/-6{x)/2A(x), with g(x) = y^WV+W) 2 , 
corresponding to the zigzag chain [12J. Within the Lan- 
dau treatment [20| , one finds that the correlation function 
of the linear chain, evaluated for a static perturbation 
at a point, decays exponentially with the length scale 
£ ~ awo/ y // S(x). In the following we focus on the situ- 
ation in which the trap frequency in the z-axis is much 
larger than the y-axis so that ifi(x) is along y and we can 
drop the a label of the field. 

Within the GL description we now assume that the 
transverse trap frequency v t undergoes a change in time 
in the interval [—tq,tq], sweeping through the mechan- 
ical instability from the linear to the zigzag chain, such 



that v 2 = v [ t c> (0) 2 -S ^ and u\ c > '(O) 2 > 5 > 0. In this 
parameter regime, we can use the time-dependent param- 
eter 5(x,t) = v 2 {t) — v[ c \x) 2 inside the GL equation. We 
also assume that the chain is in contact with a thermal 
reservoir at low temperature T, which is warranted by 
laser cooling the chain motion. More specifically, we as- 
sume that some ions of the chain are Doppler cooled. In 
the Lamb-Dicke regime, where the mechanical effects of 
atom-photon interactions can be treated in perturbation 
theory, the energy distribution of the crystal modes obeys 
a Fokker-Planck equation [2l|. The equation of motion 
for the field can be then written as 

d 2 ip - h{x) 2 d 2 ip + f]d t ^ + 6(x, t)ip + 2A{x)ip 3 = e{t) (2) 
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where the scalar e(t) is the Langevin force, describing the 
diffusion due to laser cooling, such that its moments ful- 
fill the relations (e(t)) = 0, (e(t)e(t')) = 2rjK B T5(t - t'), 
where Kg is Boltzmann constant and T is the tempera- 
ture of Doppler cooling (22[. In deriving Eq. ([5]), we have 
neglected axial distortions of the charge density due to 
the value of v t . 



We now estimate defect formation following a quench 
in the tranverse trapping frequency [23| . The nucle- 
ation of defects in such scenario resembles the forma- 
tion of solitons in a cigar-shaped Bose-Einstein con- 
densate recently discussed by Zurek As a re- 
sult of the inhomogneous charge distribution in the sys- 
tem, the transverse frequency is quenched through the 
critical point at different times along the chain, giv- 
ing rise to a propagating front along the axis, whose 
coordinates [xp,tp) satisfy 5(xF,tp) — 0. The front 
velocity, at which the instability propagates, can be 
found by taking the ratio between the characteristic 
length of the control parameter, (d x 5(x, t)/S(x, t))~ , 
and the characteristic time scale at which it changes, 
(d t 6(x,t)/8(x,t)y 1 , giving v F ~ dfs(l'!t) ■ For the s P a " 
tial dependence of the local critical frequency v%(x) — 
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dx 



1 3 the front velocity takes the form vf ~ 
L5 ° 1 with X = x/L. 



Whenever the transition is homogeneous, vf becomes 
infinite and the standard scenario of defect formation 
of Kibble-Zurek applies: the density of defects in this 
case is simply determined by the correlation length at 
the freeze-out time scale t. Elsewhere, the sound veloc- 
ity comes into play. In order to compute it, we first 
note that the relative frequency can be written with 
reference to tF as 5(x,t) = —6o[t — tp(x)]/rQ, where 



tp(x) = tq[v[°\q) 2 — v\ C! (x)' 2 ]/ 8q. One can find the time 
scale relative to tp, t, at which the dynamics stop be- 
ing adiabatic by equating the time scale 5/8 to the re- 
laxation time t. Two regimes can be identified, which 
refer to the relation between the damping ratio and the 
value of 8 at i. The so-called overdamped regime [5[ 

corresponds to the situation r\ 

case, t = (titq/8o) 1/2 , 

relation length £ x 



I8(p,t). In this 
which sets the freezed-out cor- 
aujQ(-q8o/T Q y 1 / i . 



aaj /y\8(x,t)\ 
Then, the characteristic velocity of a perturbation be- 
comes v x = £ x /t x — aLOo^o/rfTQ) 1 / 2 . The condition 
for kinks formation reads ^r- — A t^j(1 — X 2 )~ 2 > I 
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One can estimate the 



with A = - {c) \ 2 , 

effective size of the chain 2X* where the homogeneous 
KZM applies, by setting ^ = 1, and assuming I, < 1 

whence it follows that \XJ ~ A - The density of kinks 



obeys then the relation 



2\X* 
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(3) 



Note that this leads to a stronger dependence on tq than 
in the homogeneous case, where defects can nucleate all 

over the system, and d Q - I" 1 = ±^ . By 

contrast, in the underdamped regime (tq -c Sq/tj 3 ), the 
relaxation time is independent of the dissipation and di- 
verges as t = 1/ y/\8(x, t)\, which leads to the freeze-out 
time scale t = (tq/Sq) 1 / 3 . At this time scale, the cor- 
relation length reads £ x = acjo(rQ/(5o) 1 ^ 3 leading to a 
uniform sound velocity v x = £,x/t z = aujQ. The causal- 
ity argument implies that -P 1 = _4„-p^(I — X 2 )~ 2 > I, 

in terms of the parameter A u = ,^,_ L „ — • For the 



purpose of deriving a scaling of the density of defects, 
we assume that formation of kinks arises only in a small 
central region X* < 1 of approximate size 2|X*| ~ A u , 
so that 
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L 



3i/ t (c) (0) 2 a 2 wg 
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which should be compared with the density of defects 

£ , , . / s \V3 
in the homogeneous case d u ~ t = -— — 

We shall refer to the mechanism above as the inhomo- 
geneous KZM (IKZM), whose main prediction is the 
scaling in Eqs. @ and in dramatic contrast with 
their homogeneous counterparts. To test the IKZM 
we next consider the dynamics of the structural phase 
transition for ions in a linear Paul trap, whose trans- 
verse square-frequency is driven through a linear quench 
as above. We consider that the phase transition is 
well described by a Langevin dynamics following [H, [2l| . 
The ions obey the set of coupled stochastic equations 
mr s: + d rz V({r,},t) + mqr, + e{t) — 0, (i — 1,...,N) 
where = (xi,yi), s(t) is Langevin force (whose ampli- 
tude we take e = 0.05 mloi/ 2 with Iq = Q 2 /mv 2 ), and r] is 
the damping constant. At t — the chain is in the ground 
state, with all the ions at the equilibrium position satis- 
fying d Ti V({ri}) = 0. The value of <5(t) is chosen such 
that the equilibrium configuration is linear, but close to 
the critical frequency below which the ground state be- 
comes doubly degenerated. The system is then driven 
through the transition. Typical defects are shown in Fig. 
[T]d, and come in two varieties, Z2 kinks (of topological 
charge a = +1) and anti-kinks (a = —I). These defects 
resemble the nonmassive kinks of the Frenkel-Kontorova 
model with a transversal degree of freedom which can be 
described by an effective </> 4 theory for the translational 
displacement [3] • When defects appear near the edges of 
the chain, they might be lost in the linear part. To min- 
imise defects losses at the edges of the chain (see below) , 
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FIG. 2: Density of defects for a harmonically trapped ion 
chain as a function of the inverse of the sweeping rate (a) in 
the overdamped regime (rj = 100^) , where the slope in the fit 
is 0.995 with regression coefficient 1. (b) in the underdamped 
regime (rj = 10^), where the slope in the fit is 1.427 with 
regression coefficient 0.994. The defects are only considered 
in the central Nc = 30 ions, in order to minimize defect losses 
(N = 50, 2000 realizations). 

the density of defects d (number of defects over the total 
number of ions) is computed once the average absolute 
transverse displacement of the ions (y) = J2iec \Vi\/Nc 
approaches 90% of that of the ground state in the final 
trap. C denotes the set of Nc central ions which would 
reach the zigzag structure in an adiabatic transition, and 
where the formation of defects is studied. The average 
density of defects d over different realizations such as the 
one in Fig. |Tfb) is computed for different values of tq, 
and a least-squares fit to the list of data is used to extract 
the exponent governing the scaling. Numerical simula- 
tions in Fig. [2] are in good agreement with the IKZM 
scaling derived in Eqs. ([3]) and Finite-size effects of 
the chain lead to a saturation of the density of defects 
and deviations from IKZM. Further, the applicability of 
IKZM is restricted by the following effects: a) Axial and 
transverse modes are coupled since the ions shift in the 
axial direction towards the center of the trap as the struc- 
tural phase transition takes place, b) The amplitude of 
the transverse displacement of the ions increases in the 
center of the trap, making the effective Peierls-Nabarro 
potential seen by a kink [2J] the convolution of a peri- 
odic potential with an inverted bell-shaped function, and 
leading to transport of defects and losses near the edges 
of the trap. Defect transport remains even if the longitu- 
dinal degrees of freedom of the ions are frozen on a lattice 
due to the transverse motion, and even when the trapping 
potential makes the inter-ion spacing homogeneous due 
to a local correction to the transverse critical frequency 
in the finite system, c) Defects with the same topological 
charge repel each other and attract otherwise. Scatter- 
ing between kinks and anti-kinks can occur leading to 
their annihilation, a process particularly relevant in the 
underdamped regime which leads to deviations from the 
IKZM. 

In conclusion, we have proposed an ion crystal as 
a test-bed for the formation of structural defects gov- 
erned by the KZM. This system is far more amenable 
to experimental verification and control than other sys- 
tems with realistic possibilities to enter the quantum 



regime. Though the paradigmatic result for an homo- 
geneous second-order phase transition can be studied us- 
ing a ring trap, the inhomogeneities in a linear Paul trap 
make an ion crystal a natural system to study the IKZM 
mechanism where the scaling of the number of topolog- 
ical defects as a function of the quenching rates is dra- 
matically altered, as we have shown analytically and con- 
firmed by numerical simulations. 
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